\(h_{0}(t)\) er baseline hazard ved tiden \(t\) for en referenceperson (en person med \(x=0\)). Baseline hazard specificeres ikke på forhånd, hvilket gør modellen semi-parametrisk. Bemærk, at der ikke er et separat intercept, fordi denne information er “absorberet” i \(h0(t)\). (Se også slides fra forrige lektion).
Cox PH: fixed covariater
Hazarden for den \(i\)’ende person til tiden \(t\), er et produkt af to led i:
\(h_{0}(t)\) er baseline hazard ved tiden \(t\) for en referenceperson (en person med \(x=0\)). Baseline hazard specificeres ikke på forhånd, hvilket gør modellen semi-parametrisk. Bemærk, at der ikke er et separat intercept, fordi denne information er “absorberet” i \(h_{0}(t)\). (Se også slides fra forrige lektion).
\(e^{\beta_{k} x_{k}}\) angiver effekten af de \(k\) covariater. Effekten er konstant over tid og multiplicerer (på) baseline hazard for at give den samlede risiko for den enkelte person med specifikke værdier for covariaterne.
Cox PH: forudsætninger
Forudsætning for Cox ’s regressionsmodel er proportionalitet mellem hazard-raterne (den proportionale hazard model).
For to grupper gælder det at:
\[
h_{i}(t) = \theta h_{j}(t)
\]
Derfor: ved sammenligning af hazard-raterne for to personer (\(i\) og \(j\)) på sammentidspunkt gælder det for hazard-ration at:
\[
\hat{HR} =
\frac{h_{i}(t)}{h_{j}(t)} = \theta
\] hvor \(\theta\) er en konstant.
Cox PH: Hazard-ration
Hazard-ratioen sammenligner risikoen for en begivenhed for en person uden for referencegruppen (\(i\)) med en person i referencegruppen (\(j\)), udtrykt som:
Bemærk at \(h_{0}(t)\) ophæves, så HR kun afhænger af \(\text{exp}(\beta)\). Dette svarer konceptuelt til, hvordan Odds Ratio (OR) fungerer i logistisk regression.
Hvis HR er \(=1\) har personen uden for referencegruppen samme risiko som personen i referencegruppen.
Hvis HR er \(<1\) er risikoen lavere for personen uden for referencegruppen.
Hvis HR er \(>1\) er isikoen højere for personen uden for referencegruppen.
Da HR er et forhold (en ratio), udtrykker den odds—ikke direkte sandsynligheder.
Cox PH: Hazard-ration og proportionalitet
I første figur ses en situation, hvor den estimerede hazard-ratio er konstant over tid, altså \(\hat{HR} = \frac{h_{i}(t)}{h_{j}(t)} = \theta\). I det andet plot varierer den estimerede hazard-ratio ikke-proportionelt over tid, så \(\hat{HR} = \frac{h_{i}(t)}{h_{j}(t)} \neq \theta\).
Den centrale pointe her er at den absolutte hazard-rate \(h(t)\) kan ændre sig over tid, fordi Cox-modellen er semi-parametrisk og kan tilpasses data med forskellige sandsynlighedstæthedsfunktioner (pdf). Det afgørende er, at forholdet mellem hazardene for de to grupper—hazard-ratioen—skal være proportional. Det betyder, at selvom den samlede risiko (hazardraten) kan variere over tid, skal variationen være ens for begge grupper, så sammenligningen mellem dem (HR) forbliver konstant.
klik for at se koden
library(tidyverse)# Tid:t <-seq(0, 10, by =0.1)# beta-værdien: effekten af covariaten, k:beta <-0.5# dataframe:df <-tibble(tid = t,h0 =0.1+0.05*sin(t), # Baseline hazardh1 = h0 *exp(beta), # Hazard for gruppe iHR = h1 / h0 # Hazard-ratio mellem i og j)# Plot:p1 <-ggplot(df, aes(x = tid)) +geom_line(aes(y = h0, color ="j"), size =1.2) +geom_line(aes(y = h1, color ="i"), size =1.2) +labs(x ="Time", y ="Hazard", title ="Proportionale Hazardfunktioner", color ="Gruppe") +theme_minimal()p1
klik for at se koden
library(tidyverse)# Tid:t <-seq(0, 10, by =0.1)# Baseline hazard-funktion for referencegruppen, j:h0 <-0.1+0.05*sin(t)# Beta (effekten) + "gamma" (en tidsafhængig effekt der varierer over tid):beta <-0.5gamma <-0.2# Hazard-funktion for gruppe i med tidsafhængig ikke-proportional effekt:h2 <- h0 *exp(beta + gamma * t)# Hazard-ratioen:HR <- h2 / h0 # Dette svarer til exp(beta + gamma*t) ... altså, jeg har lavet en teoretisk funktion, der er ikke-proportionel over tid:# Dataframe:df <-tibble(tid = t,h0 = h0,h2 = h2,HR = HR)# Plot:p1 <-ggplot(df, aes(x = tid)) +geom_line(aes(y = h0, color ="j"), size =1.2) +geom_line(aes(y = h2, color ="i"), size =1.2) +labs(x ="Tid", y ="Hazard", title ="Ikke-Proportionale Hazardfunktioner", color ="Gruppe") +theme_minimal()p1
Cox PH: Hazard-ration og proportionalitet
I første figur ses en situation, hvor den estimerede hazard-ratio er konstant over tid, altså \(\hat{HR} = \frac{h_{i}(t)}{h_{j}(t)} = \theta\). I det andet plot varierer den estimerede hazard-ratio ikke-proportionelt over tid, så \(\hat{HR} = \frac{h_{i}(t)}{h_{j}(t)} \neq \theta\).
Bemærk også at proportionalitet betyder, at forholdet (HR) mellem de to hazardfunktioner forbliver konstant over tid (i.e., den ekstra risiko, der tilføjes af covariat-effekten, \(\exp(\beta)\), i forhold til baseline), men det betyder ikke, at afstanden (den absolutte forskel) mellem dem er konstant. Eftersom \(h_{0}(t)\) kan variere over tid, kan den absolutte forskel også variere. Derfor: selvom \(\beta\) bliver eksponentieret og dermed giver et konstant forhold mellem grupper, kan afstanden mellem linjerne (funktionerne) være større, når den underliggende baseline hazard er højere.
klik for at se koden
library(tidyverse)# Tid:t <-seq(0, 10, by =0.1)# beta-værdien: effekten af covariaten, k:beta <-0.5# dataframe:df <-tibble(tid = t,h0 =0.1+0.05*sin(t), # Baseline hazardh1 = h0 *exp(beta), # Hazard for gruppe iHR = h1 / h0 # Hazard-ratio mellem i og j)# Plot:p1 <-ggplot(df, aes(x = tid)) +geom_line(aes(y = h0, color ="j"), size =1.2) +geom_line(aes(y = h1, color ="i"), size =1.2) +labs(x ="Tid", y ="Hazard", title ="Proportionale Hazardfunktioner", color ="Gruppe") +theme_minimal()p1
klik for at se koden
# Plot:p2 <-ggplot(df, aes(x = tid, y = HR)) +geom_line(color ="green", size =1.2) +labs(x ="Time", y ="Hazard Ratio", title ="Hazard Ratio (Konstant over tid)") +theme_minimal()p2
klik for at se koden
library(tidyverse)# Tid:t <-seq(0, 10, by =0.1)# Beta (effekten) + "gamma" (en tidsafhængig effekt der varierer over tid):beta <-0.5gamma <-0.2# Dataframe:df <-tibble(tid = t,h0 =0.1+0.05*sin(t),h2 = h0 *exp(beta + gamma * t), # Hazard-funktion for gruppe i med tidsafhængig ikke-proportional effektHR = h2 / h0 # Dette svarer til exp(beta + gamma*t) ... altså, jeg har lavet en teoretisk funktion, der er ikke-proportionel over tid)# Plot:p1 <-ggplot(df, aes(x = tid)) +geom_line(aes(y = h0, color ="j"), size =1.2) +geom_line(aes(y = h2, color ="i"), size =1.2) +labs(x ="Tid", y ="Hazard", title ="Ikke-Proportionale Hazardfunktioner", color ="Gruppe") +theme_minimal()p1
klik for at se koden
# Plot:p2 <-ggplot(df, aes(x = tid, y = HR)) +geom_line(color ="green", size =1.2) +labs(x ="Tid", y ="Hazard Ratio", title ="Hazard Ratio (Ikke-konstant over tid)") +theme_minimal()p2
Tied data
Tied data opstår ofte i overlevelsesanalyser, fordi vi typisk måler tid i diskrete enheder. Det betyder, at flere hændelser kan ske i samme tidsinterval. Denne situation bryder med antagelsen om, at alle observationer er uafhængige, og kan derfor føre til bias i estimater og standardfejl, som påvirker validiteten af vores statistiske konklusioner.
Vi kan løse dette problem ved hjælp af tre forskellige standard-metoder:
Breslow: Velegnet til “små” datasæt med få ties.
Efron: Klarer bedre mange ties.
Exact: Den mest nøjagtige metode, men også den mest beregningstunge (computerkrævende).
Tied data: Efron
Når flere personer oplever hændelsen i samme tidsinterval, kaldes det for tied data. Efron-metoden håndterer dette ved at tage højde for alle mulige rækkefølger (rank orderings) af de hændelser, der sker på samme tid. På den måde justeres, hvor meget de “tied” hændelser bidrager til den samlede likelihood-funktion når vi estimerer vores model.
Vi udtrykke med denne metode udtrykke hazardfunktion ved en distinkt tid \(t_{(i)}\) med tied events som:
\(h_{a}(t_{(i)})\) er den kontrollerede/justerede hazard funktion ved \(t_{(i)}\),
\(t_{(i)}\) er den distinkte tid, hvor hændelser opstå,
\(j\) er indekset for individer som oplevede begivenheden ved \(t_{(i)}\),
\(d_{i}\) er antallet af individer der oplevede begivenheden ved tid \(t_{(i)}\),
\(R(t_{(i)})\) er risikosættet før \(t_{(i)}\),
\(k\) er indekset for individer \(i\)risikosættet, \(R(t_{(i)})\), før \(t_{(i)}\),
\(\delta_{ij}\) er en dummy, der indikerer hvis individ \(j\) oplevede begivenheden ved \(t_{(i)}\) (1) eller ej (0),
\(X_{j}\) er en vector af covariater for individ \(j\),
\(\beta\) er en vector af coeficcienter fra Cox PH modellen.
Tied data: Efron
Metoden tager højde for alle mulige rækkefølger af disse forbundne begivenheder inden for et givent tidsinterval. For eksempel:
3 personer, \(A, B, C\), oplevede begivenheden ved \(t=x\).
Da vi arbejder med diskrete tidsenheder er der en række af mulige “korrekte” rækkefølger som \(A\), \(B\) og \(C\) oplevede begivenheden:
\(A\) først, så \(B\), så \(C\)
\(B\), så \(C\), så \(A\),
osv. …
På den måde “glatter”/smoother vi estimatet af hazardfunktionen og kontrollerer for den usikkerhed, der opstår, når begivenheder registreres som værende samtidige og observationerne potentielt er afhængige.
Tied data: Efron
Det vil sige, at vi overvejer alle mulige måder, de individer, der oplever begivenheden samtidigt, kan rangeres, hvilket justerer for den usikkerhed, der opstår ved tied data. Derfor får vi en mere præcis likelihood funktion, udtrykt som:
hvor produktet over \(i\) går fra 1 to \(m\). \(m\) er det samlede antal diskrete begivenhedstindspunkter.
Fordi vi kun tager de individer, \(j\), der oplevede begivenheden i betragtning og ikke specificerer baseline hazard funktionen er det formelt en partial likelihood funktion.
Tied data: Exact
På mange måder ligner Exact metode Efron metoder. Den grundlæggende forskel er præcisionen i estimater (og den direkte relaterede computerkraft det kræver at være mere præcis). Med andre ord er \(h(t)\) ikke bare “kontrollere”/adjusted som med Efron metoden, men præcis.
Hvor Efron kontrollerede for mulige permutationer af forbundne begivenheder indenfor et sæt af forbundne begivenheder, tager “exact” metode alle mulige permutationer i dataen (og ikke bare sættet) i betragtning og er derfor ikke afhængig af et “smoothed” estimat men giver et præcist estimat.
Cox PH med tidsvariende variable
For at udvide Cox modellen med tidsvarierende covariate er der kun ganske få ændringer i formlen I er ved at være bekendt med:
Den eneste forskel er at vi tilføjer \(t\) for de tidsvarierende covariater.
Her er \(x_{2}(t)\) en tidsvarierende variabel, mens de øvrige \(x\)-variable er faste, og skal læses som at hazarden ved tid \(t\) er bestemt af tidsinvariante karakteristika, \(x_{1}\), og \(x_{2}\)ved tid\(t\).
Det betyder, at hazard-ratioen i denne situation bliver tidsafhængig, fordi den afhænger af, hvordan den tidsvarierende variabel adskiller sig mellem personerne på et givet tidspunkt.
Cox PH med tidsvariende variable: proportionalitet?
klik for at se koden
library(tidyverse)# Tidt <-seq(0, 10, by =0.1)# Beta: effekten af covariatenbeta <-0.5# Tidsvarierende covariater for de to grupper:# Gruppe i: x_i(t) = 1 + 0.1*t# Gruppe j: x_j(t) = 1 + 0.05*tx_i <-1+0.1* tx_j <-1+0.05* t# Beregn hazard for hver gruppe:h_i <- h0 *exp(beta * x_i) # Hazard for gruppe ih_j <- h0 *exp(beta * x_j) # Hazard for gruppe j# Opret en tibble med datadf <-tibble(tid = t,h0 =0.1+0.05*sin(t),`i`= h_i,`j`= h_j,HR =exp(beta * (x_i - x_j)))# Omform data til long format til ggplotdf_long <- df %>%pivot_longer(cols =c("i", "j"),names_to ="Gruppe",values_to ="Hazard")# Plot:p1 <-ggplot(df_long, aes(x = tid, y = Hazard, color = Gruppe)) +geom_line(size =1.2) +labs(title ="Hazardfunktioner over tid for to grupper",x ="Tid",y ="Hazard") +theme_minimal()p1
klik for at se koden
library(tidyverse)# Tidt <-seq(0, 10, by =0.1)# Beta: effekten af covariatenbeta <-0.5# Tidsvarierende covariater for de to grupper:# Gruppe i: x_i(t) = 1 + 0.1*t# Gruppe j: x_j(t) = 1 + 0.05*tx_i <-1+0.1* tx_j <-1+0.05* t# Beregn hazard for hver gruppe:h_i <- h0 *exp(beta * x_i) # Hazard for gruppe ih_j <- h0 *exp(beta * x_j) # Hazard for gruppe j# Opret en tibble med datadf <-tibble(tid = t,h0 =0.1+0.05*sin(t),`i`= h_i,`j`= h_j,HR =exp(beta * (x_i - x_j)))# Omform data til long format til ggplotdf_long <- df %>%pivot_longer(cols =c("i", "j"),names_to ="Gruppe",values_to ="Hazard")# Plot: p2 <-ggplot(df, aes(x = tid, y = HR)) +geom_line(color ="green", size =1.2) +labs(title ="Tidsvarierende Hazard Ratio",x ="Tid",y ="Hazard Ratio, HR(t)" ) +theme_minimal()p2
Cox PH med tidsvariende variable: Data
flowchart TD
X([<font size=6><font color=white> Afhængig <br/> variabel])
A([<font size=6><font color=black> Forklarende <br/> variabel])
-...-> B([<font size=6><font color=white>Tidsvariaende])
A
---> C([<font size=6><font color=black>Konstant])
B
--> D([<font size=6><font color=black> Variable der løber <br/> med tiden])
B
--> E([<font size=6><font color=black>En variabel med forskellige <br/> værdier til forskellige <br/> tidspunkter])
style A fill:#f45c2c,stroke:#111,stroke-width:1px
style B fill:#333,stroke:#111,stroke-width:1px
style C fill:#d8897b,stroke:#111,stroke-width:1px
style D fill:#9f6671,stroke:#111,stroke-width:1px
style E fill:#5e4f6d,stroke:#111,stroke-width:1px
Cox PH med tidsvariende variable: Data
flowchart TD
B([<font size=6><font color=white>Tidsvariaende])
--> D([<font size=6><font color=black> Variable der løber <br/> med tiden])
B
--> E([<font size=6><font color=black>En variabel med forskellige <br/> værdier til forskellige <br/> tidspunkter])
D
--- F([For eksempel alder])
D
-...- G([<font size=6><font color=black>Kan altid anvendes bestemmes <br/> ved start- og slutoplysninger])
E
--- H([For eksempel arb. eller udd. status])
E
-...- I([<font size=6><font color=black>Skal bestemmes for <br/> hvert tidspunkt])
style B fill:#d8897b,stroke:#111,stroke-width:1px
style D fill:#9f6671,stroke:#111,stroke-width:1px
style E fill:#5e4f6d,stroke:#111,stroke-width:1px
style G fill:#b7bd50,stroke:#111,stroke-width:1px
style I fill:#b7bd50,stroke:#111,stroke-width:1px
Cox PH med tidsvariende variable: Data
For at coxph() funktionen kan tage tidsvarierende informationer er det (desværre) ikke nok med periode_start og periode_slut. For hvert tidspunkt og person skal vi have start og slut information.
Konstruér subject-period file
Fiktive data
id
år
børn
alder
region18
tid
event
person_start
person_slut
begivenhed
periode_slut
1
1
2000
0
18
Nord
1
0
1
0
0
0
1
2001
0
19
Nord
2
0
0
0
0
0
1
2002
1
20
Syd
3
1
0
0
1
1
2
2
2003
0
20
Syd
3
0
1
0
0
0
2
2004
0
21
Øst
4
0
0
0
0
0
2
2005
0
22
Øst
5
0
0
0
0
0
2
2006
0
23
Nord
6
0
0
0
0
0
2
2007
0
24
Nord
7
0
0
1
0
1
3
3
1999
0
18
Nord
1
0
1
0
0
0
3
2000
0
19
Nord
2
0
0
0
0
0
3
2001
0
20
Vest
3
0
0
0
0
0
3
2002
0
21
Vest
4
0
0
0
0
0
3
2003
1
22
Nord
5
1
0
1
1
1
4
4
1993
0
23
Øst
6
0
1
0
0
0
4
1994
0
24
Øst
7
0
0
0
0
0
4
1995
0
25
Øst
8
0
0
0
0
0
4
1996
1
26
Øst
9
1
0
0
1
1
Cox PH med tidsvariende variable: Data
For at coxph() funktionen kan tage tidsvarierende informationer er det (desværre) ikke nok med periode_start og periode_slut. For hvert tidspunkt og person skal vi have start og slut information.
Konstruér subject-period file
Fiktive data
id
år
børn
alder
region18
tid
event
person_start
person_slut
begivenhed
periode_slut
start
slut
1
1
2000
0
18
Nord
1
0
1
0
0
0
1
2
1
2001
0
19
Nord
2
0
0
0
0
0
2
3
1
2002
1
20
Syd
3
1
0
0
1
1
3
4
2
2
2003
0
20
Syd
3
0
1
0
0
0
1
2
2
2004
0
21
Øst
4
0
0
0
0
0
2
3
2
2005
0
22
Øst
5
0
0
0
0
0
3
4
2
2006
0
23
Nord
6
0
0
0
0
0
4
5
2
2007
0
24
Nord
7
0
0
1
0
1
5
6
3
3
1999
0
18
Nord
1
0
1
0
0
0
1
2
3
2000
0
19
Nord
2
0
0
0
0
0
2
3
3
2001
0
20
Vest
3
0
0
0
0
0
3
4
3
2002
0
21
Vest
4
0
0
0
0
0
4
5
3
2003
1
22
Nord
5
1
0
1
1
1
5
6
4
4
1993
0
23
Øst
6
0
1
0
0
0
1
2
4
1994
0
24
Øst
7
0
0
0
0
0
2
3
4
1995
0
25
Øst
8
0
0
0
0
0
3
4
4
1996
1
26
Øst
9
1
0
0
1
1
4
5
Cox PH med tidsvariende variable: Data
For at coxph() funktionen kan tage tidsvarierende informationer er det (desværre) ikke nok med periode_start og periode_slut. For hvert tidspunkt og person skal vi have start og slut information.
library(tidyverse)library(readxl)df.rå <-read_excel("eksempel_data.xlsx")df.lagged <- df.rå %>%mutate(person_start =case_when(lag(id) == id ~0, TRUE~1),person_slut =case_when(lead(id) == id ~0, TRUE~1)) %>%group_by(id) %>%mutate(begivenhed =if_else(lag(børn) ==0& børn ==1& person_start ==0, 1, 0),periode_slut =case_when(begivenhed ==1| børn ==0& person_slut ==1~1, TRUE~0)) %>%slice(1:which.max(periode_slut ==1)) %>%mutate(start =1:n(),slut = start +1) %>%ungroup()
Cox PH med tidsvariende variable: Data
Disse to variable gør det muligt at udvide coxph() funktionen med:
Dette kalder vi også for et Counting process argument.
Her angives start og sluttiden for det enkelte interval, samt event (1=hændelse og 0=censurering). Intervallets størrelse afgøres af hvor hyppigt de tidsvarierende variable kan skifte værdi. Værdien skal være konstant inden for det enkelte interval. Denne lille ændring gør det muligt at håndtere tidsvarierende variable i cox-regressionen
Oprindeligt: - 1 % af alle individer fundet i IDAP igennem perioden 1980-2010
Reduktion: - Data indeholder personer, der er fyldt 10 år mellem 1980-1990 - Personer, der får børn før de er fyldt 18 år er taget ud af data (56 personer) - Kun personer, der er til stede i DK, det år de bliver 18 år. - Årlige observationer af disse personer
Øvelse - Første barn, fortsat
Diskutér de forskellige variable i datasættet:
Er de tidsvarierende eller ej?
Vil det give mening at behandle nogle af variablene som faste (værdien på en tidsvarierende variable ved forløbets start)
Gør Forlob-datasættet klar til analyse
Vær opmærksom på variabeltyperne (chr, num og factor) og hvordan I bruger dem i analysen
ALDER er en character variable
Sæt udgangspunktet for forløbet til at være året, hvor personerne er 18 år
Smid irrelevante observationsrække ud af datasættet (relevante obs er fra det 18. år frem til og med intervallet, hvor hændelsen sker eller observationen af personen slutter)
Vælg minimum tre variable til analysen af problemstillingen (minimum 1 tidsvarierende)
Lav indledende KM analyser for at undersøge forskelle mellem grupper.
– plot KM-analyserne
– Husk at begrænse datasættet til den sidste observation!
Lav en cox-regression, hvor tidsvarierende variable (mindst en) indgår, og fortolk resultaterne.
Diskuter outputtet. Hvad giver mest mening at tolke på, coef eller exp(coef), hvorfor?